////////////////////////////////////////////////////////////////////////////////
// Includes
////////////////////////////////////////////////////////////////////////////////

// Header
#include "TSPPlugin.hpp"

// Random number generator
#include <cstdlib>

// Vector
#include <vector>

// Sort algorithm
#include <algorithm>

// String manipulation
#include <string>
#include <sstream>

// I/O stream
#include <fstream>
#include <iostream>

// Math
#include <cmath>

// Conversion functions
#include "Convert.hpp"

// Tokeniser
#include "Tokeniser.hpp"

// Exceptions
#include "IOFileException.hpp"
#include "SizeException.hpp"
#include "AssertionException.hpp"
#include "InterruptionException.hpp"

// Namespaces
using namespace std;
using namespace Toolkit;
using namespace Types;
using namespace Exceptions;
using namespace SC;

////////////////////////////////////////////////////////////////////////////////
// Plugin Data
////////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////
////////////// Input data ///////////////
/////////////////////////////////////////

// Calling SC Program's name
static string program_name;

// Number of bits per city (max 2^BPC cities can be encoded)
static unsigned short nb_bpc;

// Distances matrix
static vector< vector<float> > distances;

// Number of crossover points
static unsigned short nb_cross_pts;

// Number of mutation points
static unsigned short nb_mutations;

// Maximum number of computations
static unsigned int nb_comp_max;

/////////////////////////////////////////
////////// Computation state ////////////
/////////////////////////////////////////

// Number or computations
static unsigned int nb_computations;

// Best distance
static unsigned int best_dist;

// Best path
static vector<unsigned int> best_path;

/////////////////////////////////////////
////////// Evolution operators //////////
/////////////////////////////////////////

// Selection type
enum SelType { KEEP_BEST, FIT_PROP };

// Operator types
enum OpType { PMX = 0, MOX, SWP, MOV, REV };

// Number of operators
#define NB_OPS 5

// Associated string names
static const string op_names[NB_OPS] = { "PMX", "MOX", "SWP", "MOV", "REV" };

// Associated operators codes (used in the SC file)
#define PMX_CODE	"000"
#define MOX_CODE	"001"
#define SWP_CODE	"010"
#define MOV_CODE	"011"
#define REV_CODE	"100"
static const string op_codes[NB_OPS] = { PMX_CODE, MOX_CODE, SWP_CODE, MOV_CODE, REV_CODE };

/////////////////////////////////////////
///////// Computation recording /////////
/////////////////////////////////////////

// Record progression
static vector< pair<unsigned int,unsigned int> > record;

// Record operations fitness
static vector< pair<float,float> > mem[NB_OPS];
static float	op_fit[NB_OPS] = { 1.f, 1.f, 1.f, 1.f, 1.f },
				op_fit_min = 1.f;

////////////////////////////////////////////////////////////////////////////////
// Tool functions
////////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////
//////////// Data functions /////////////
/////////////////////////////////////////

// Calulate the total distance of a path
static const unsigned int getTotalDist(const vector<unsigned int> & path) {
	if (path.size() == 0) return 0;
	float total_dist = distances[0][path[0]+1];
	for (unsigned int i=0; i<path.size()-1; ++i)
		total_dist += distances[path[i]+1][path[i+1]+1];
	total_dist += distances[path[path.size()-1]+1][0];
	++nb_computations;
	return (unsigned int)total_dist;
}

// Encode a path into a binary string
static const string getStringPath(
	const vector<unsigned int> & path,
	const unsigned int length
) {
	string str_path;
	// Convert the path
	for (unsigned int i=0; i<path.size(); ++i)
		str_path += Convert::UIntToBinaryString(path[i],nb_bpc);
	// Test that it is not longer than the required length
	if (str_path.size() > length) {
		stringstream ss;
		ss << "Path string length is too long (" << str_path.size() << ") to match the required length of " << length;
		throw SizeException(ss.str());
	}
	// Fill in with '?'
	while (str_path.size() < length) str_path += '?';
	// Return path string
	return str_path;
}

// Decode a path from a binary string
static const vector<unsigned int> getPathFromString(const string & str_path) {
	vector<unsigned int> path;
	for (unsigned int i=0; i<nb_bpc*(distances.size()-1); i+=nb_bpc)
		path.push_back(Convert::BinaryStringToUInt(str_path.substr(i,nb_bpc)));
	return path;
}

/////////////////////////////////////////
////////////// Operators/ ///////////////
/////////////////////////////////////////

// Mutate a path with "Swap"
static const vector<unsigned int> swap(const vector<unsigned int> & path) {
	vector<unsigned int> new_path = path;
	for (unsigned int m=0; m<nb_mutations; ++m) {
		unsigned int	idx = rand() % new_path.size(),
						tmp = new_path[idx];
		new_path[idx] = new_path[ (idx+1) % new_path.size() ];
		new_path[ (idx+1) % new_path.size() ] = tmp;
	}
	return new_path;
}

// Mutate a path with "Reverse" (2opt)
static const vector<unsigned int> reverse(const vector<unsigned int> & path) {
	vector<unsigned int> new_path = path;
	for (unsigned int m=0; m<nb_mutations; ++m) {
		unsigned int idx1 = rand() % new_path.size(), idx2 = idx1, tmp;
		while (idx2 == idx1) idx2 = rand() % new_path.size();
		if (idx2 < idx1) {
			tmp = idx1;
			idx1 = idx2;
			idx2 = tmp;
		}
		for (unsigned int j=0; j<(idx2-idx1+1)/2; ++j) {
			tmp = new_path[idx1+j];
			new_path[idx1+j] = new_path[idx2-j];
			new_path[idx2-j] = tmp;
		}
	}
	return new_path;
}

// Mutate a path with "Move"
static const vector<unsigned int> move(const vector<unsigned int> & path) {
	vector<unsigned int> new_path = path;
	for (unsigned int m=0; m<nb_mutations; ++m) {
		unsigned int idx1 = rand()%new_path.size(), idx2 = idx1, tmp;
		while (idx2 == idx1) idx2 = rand() % new_path.size();
		tmp = new_path[idx1];
		new_path[idx1] = new_path[idx2];
		new_path[idx2] = tmp;
	}
	return new_path;
}

// Partially Mapped Crossover (PMX)
static const vector<unsigned int> pmx(
	const vector<unsigned int> & keep,
	const vector<unsigned int> & replace
) {
	// Get random cross points
	vector<unsigned int> cp;
	unsigned int idx;
	for (unsigned int i=0; i<nb_cross_pts; ++i) {
		do idx = rand()%keep.size();
		while (find(cp.begin(),cp.end(),idx) != cp.end());
		cp.push_back(idx);
	}
	if ((cp.size()&1) == 1) cp.push_back(keep.size());
	sort(cp.begin(),cp.end());
	// Create new paths from parents
	vector<unsigned int> new_path = replace;
	// Copy inside the intervals from the parent not used for initialisation
	vector<unsigned int>::iterator itor;
	for (unsigned int k=0; k<cp.size()-1; k+=2) {
		for (unsigned int p=cp[k]; p<cp[k+1]; ++p) {
			if ((itor = find(new_path.begin(),new_path.end(),keep[p])) == new_path.end()) {
				stringstream ss;
				ss << "Cannot read symbol " << keep[p] << " from new_path";
				throw Exception(ss.str());
			}
			*itor = new_path[p];
			new_path[p] = keep[p];
		}
	}
	// Return the new path
	return new_path;
}

// Multi-points Order Crossover (MpOX)
static const vector<unsigned int> mox(
	const vector<unsigned int> & keep,
	const vector<unsigned int> & replace
) {
	// Get random cross points
	vector<unsigned int> cp;
	unsigned int idx;
	for (unsigned int i=0; i<nb_cross_pts; ++i) {
		do idx = rand()%keep.size();
		while (find(cp.begin(),cp.end(),idx) != cp.end());
		cp.push_back(idx);
	}
	if ((cp.size()&1) == 1) cp.push_back(keep.size());
	sort(cp.begin(),cp.end());
	// Initialise paths
	vector<unsigned int> new_path(keep.size(),0xFFFFFFFF);
	// Copy inside the intervals from one parent
	for (unsigned int k=0; k<cp.size()-1; k+=2) {
		for (unsigned int p=cp[k]; p<cp[k+1]; ++p)
			new_path[p] = keep[p];
	}
	// Complete what is outside the intervals from the other parent
	vector<unsigned int>::iterator itor;
	unsigned int idx1 = 0, lb, ub;
	for (int k=-1; k<(int)cp.size(); k+=2) {
		if (k==-1) {
			lb = 0;
			ub = cp[0];
		}
		else if (k==((int)cp.size())-1) {
			lb = cp[cp.size()-1];
			ub = replace.size();
		}
		else {
			lb = cp[k];
			ub = cp[k+1];
		}
		for (unsigned int p=lb; p<ub; ++p) {
			while (find(new_path.begin(),new_path.end(),replace[idx1]) != new_path.end()) ++idx1;
			new_path[p] = replace[idx1];
		}
	}
	// Return the new path
	return new_path;
}

/////////////////////////////////////////
//////// Selection and evolution ////////
/////////////////////////////////////////

static const bool selectAndEvolve (
	SCSystem & sys1,
	const SCSystem & context,
	SCSystem & sys2,
	const SelType sel_type,
	const OpType op_type
) {
	// Get the path and the total distance of the systems
	vector<unsigned int>	path1 = getPathFromString(sys1.GetSchema1()),
							path2 = getPathFromString(sys2.GetSchema1());
	if (path1.size() != distances.size()-1) {
		stringstream ss;
		ss << "Path1 should contain " << (distances.size()-1) << " cities but contains " << path1.size();
		throw SizeException(ss.str());
	}
	if (path2.size() != distances.size()-1) {
		stringstream ss;
		ss << "Path2 should contain " << (distances.size()-1) << " cities but contains " << path2.size();
		throw SizeException(ss.str());
	}
	unsigned int	dist1 = Convert::BinaryStringToUInt(sys1.GetSchema2()),
					dist2 = Convert::BinaryStringToUInt(sys2.GetSchema2()),
					r_dist;
	// Selection
	vector<unsigned int> *keep, *replace;
	if (sel_type == KEEP_BEST) {
		// Keep the best path and replace the other one
		if ( (dist1 < dist2) || ( (dist1 == dist2) && (rand()&1) ) ) {
			keep	= &path1;
			replace	= &path2;	r_dist = dist2;
		}
		else {
			keep	= &path2;
			replace	= &path1;	r_dist = dist1;
		}
	}
	else { // Fitness Proportional
		// Probability that sys1 is kept
		if ( (((float)rand())/((float)RAND_MAX)) <= (((float)dist2)/((float)(dist1+dist2))) ) {
			keep	= &path1;
			replace	= &path2;	r_dist = dist2;
		}
		// Otherwise keep sys2
		else {
			keep	= &path2;
			replace	= &path1;	r_dist = dist1;
		}
	}
	// Evolution (call the appropriate operator)
	if (op_type == PMX)
		*replace = pmx(*keep,*replace);
	else if (op_type == MOX)
		*replace = mox(*keep,*replace);
	else if (op_type == SWP)
		*replace = swap(*keep);
	else if (op_type == MOV)
		*replace = move(*keep);
	else
		*replace = reverse(*keep);
	// Compute the new distance
	unsigned int tmp_dist = getTotalDist(*replace);
	if (tmp_dist == 0) throw AssertionException("Path distance cannot be zero");
	// Penalise copies (giving 0.9 instead of 1.0) to maintain diversity
	if ((tmp_dist == r_dist) && (path1 == path2))
		mem[op_type].push_back(pair<float,float>((float)nb_computations,0.9f));
	else {
		mem[op_type].push_back(pair<float,float>((float)nb_computations,((float)r_dist)/((float)tmp_dist)));
		r_dist = tmp_dist;
	}
	// Store data into the appropriate system
	if (replace == &path1) {
		sys1.ReplaceSchema1(getStringPath(*replace,sys1.GetSchema1().size()));
		sys1.ReplaceSchema2(Convert::UIntToBinaryString(r_dist,sys1.GetSchema2().size()));
	}
	else {
		sys2.ReplaceSchema1(getStringPath(*replace,sys2.GetSchema1().size()));
		sys2.ReplaceSchema2(Convert::UIntToBinaryString(r_dist,sys2.GetSchema2().size()));
	}
	// Keep note of the best path
	if (r_dist < best_dist) {
		best_dist = r_dist;
		best_path = *replace;
		cout << "New best path discovered with " << sel_type << " " << op_type << endl;
		cout << "\t1-";
		for (unsigned int i=0; i<best_path.size()-1; ++i)
			cout << (best_path[i]+2) << "-";
		cout << (best_path[best_path.size()-1]+2) << endl;
		cout << "\tNb fitness computations : " << nb_computations << endl;
		cout << "\tBest solution's length  : " << best_dist << endl;
		record.push_back(pair<unsigned int,unsigned int>(nb_computations,best_dist));
	}
	// Return OK
	if ((nb_comp_max == 0) || (nb_computations < nb_comp_max)) return true;
	else {
		stringstream ss;
		ss << "Maximum number of computations (" << nb_comp_max << ") has been reached - Computation Complete";
		throw InterruptionException(ss.str());
	}
}

/////////////////////////////////////////
//////////// Self-Adaptation ////////////
/////////////////////////////////////////

// Varying window size
#define ADAPT_WINDOW_SIZE(N) ((int)(1000.f*(1.f-exp(-0.0001f*N))))
//((int)(10*sqrt(nb_computations)))

// Minimum number of samples required to compute an average
#define ADAPT_MIN_SAMPLE_SIZE 5

// Return the fitness value of the given operator
static const float getOpFit(const string & op) {
	OpType type;
	// Get type
	if (op == PMX_CODE)			type = PMX;
	else if (op == MOX_CODE)	type = MOX;
	else if (op == SWP_CODE)	type = SWP;
	else if (op == MOV_CODE)	type = MOV;
	else if (op == REV_CODE)	type = REV;
	else throw AssertionException("Operator identification error");
	int idx = ((int)mem[type].size()) - 1;
	int nb_comps = ADAPT_WINDOW_SIZE(nb_computations);
	while ((idx >= 0) && (((int)mem[type][idx].first) >= ((int)nb_computations) - nb_comps)) --idx;
	if (((int)mem[type].size())-1-idx >= ADAPT_MIN_SAMPLE_SIZE) {
		float ave = 0.f;
		for (int i=idx+1; i<(int)mem[type].size(); ++i)
			ave += mem[type][i].second;
		ave /= ((float)(((int)mem[type].size())-1-idx));
		op_fit[type] = ave;
		if (ave < op_fit_min) op_fit_min = ave;
		return ave;
	}
	else
		return -1.f;
}

// Wheel bias (no bias seems to be better)
#define WHEEL_BIAS(FIT) (FIT) // - op_fit_min/2.f)

// Wheel method to pick up an operator
static const unsigned short op_wheel() {
	float tmp = 0.f, wheel;
	for (unsigned short i=0; i<NB_OPS; ++i) tmp += WHEEL_BIAS(op_fit[i]);
	wheel = ((float)rand()*tmp)/((float)RAND_MAX);
	short idx = 0;
	tmp = WHEEL_BIAS(op_fit[idx]);
	while ((idx<NB_OPS-1) && (tmp<wheel)) {
		++idx;
		tmp += WHEEL_BIAS(op_fit[idx]);
	}
	return idx;
}

////////////////////////////////////////////////////////////////////////////////
// Plugin init and exit
////////////////////////////////////////////////////////////////////////////////

void InitPlugin (const string & prog_name) {
	// Reset data
	program_name = prog_name;
	best_dist = 0xFFFFFFFF;
	best_path.clear();
	distances.clear();
	nb_computations = 0;
	nb_cross_pts = 2;
	nb_mutations = 1;
	record.clear();
	// Input file stream
	ifstream data;
	// ----- Data File -----
	// Load the data file
	string tsp_data = "tsp_data.dat";
	data.open(tsp_data.c_str());
	// Test the file
	if (!data.is_open()) throw IOFileException("Cannot read from " + tsp_data);
	// Read data
	try {	
		// Read the number of cities
		unsigned int nb_cities;
		if (data.eof()) throw IOFileException(tsp_data + " : Cannot read the number of cities");
		data >> nb_cities;
		// Read the distances matrix
		distances.resize(nb_cities);
		for (unsigned int l=0; l<nb_cities; ++l)
			distances[l].resize(nb_cities,0);
		for (unsigned int l=0; l<nb_cities; ++l) {
			for (unsigned int c=0; c<nb_cities; ++c) {
				if (data.eof()) {
					stringstream ss;
					ss << tsp_data + " : Cannot read the distance between cities " << l << " and " << c;
					throw IOFileException(ss.str());
				}
				data >> distances[l][c];
			}
		}
		cout << "=====================================" << endl;
		cout << "=             TSP Data              =" << endl;
		cout << "=====================================" << endl;
		cout << "Number of cities : " << nb_cities << endl;
		cout << "=====================================" << endl;
		// Close file
		data.close();
	}
	catch (exception &) {
		data.close();
		throw;
	}
	// ----- Setup File -----
	data.clear();
	// Load the data file
	string tsp_setup = "tsp_setup.ini";
	data.open(tsp_setup.c_str());
	// Test the file
	if (!data.is_open()) throw IOFileException("Cannot read from " + tsp_setup);
	// Read data
	try {
		Tokeniser tokeniser(data);
		tokeniser.AddSkipSymbol(" ");
		tokeniser.AddSkipSymbol("\t");
		tokeniser.AddSkipSymbol("\n");
		tokeniser.AddCommentLineSymbol("//");
		bool bpc_ok = false, co_ok = false, mut_ok = false, nbc_ok = false;;
		while (!bpc_ok || !co_ok || !mut_ok || !nbc_ok || !data.eof()) {
			string token = tokeniser.GetNextToken(false);
			if (token == "bits_per_city") {
				if (bpc_ok) throw IOFileException("The number of bits per city is defined more than once");
				token = tokeniser.GetNextToken(true);
				nb_bpc = (unsigned short)Convert::StringToInt(token);
				if (nb_bpc == 0) throw IOFileException("The number of bits per city must be positive");
				bpc_ok = true;
			}
			else if (token == "crossover_points") {
				if (co_ok) throw IOFileException("The number of crossover points is defined more than once");
				token = tokeniser.GetNextToken(true);
				nb_cross_pts = (unsigned short)Convert::StringToInt(token);
				if (nb_cross_pts == 0) throw IOFileException("The number of crossover points must be positive");
				co_ok = true;
			}
			else if (token == "mutations") {
				if (mut_ok) throw IOFileException("The number of mutation points is defined more than once");
				token = tokeniser.GetNextToken(true);
				nb_mutations = (unsigned short)Convert::StringToInt(token);
				if (nb_mutations == 0) throw IOFileException("The number of mutation points must be positive");
				mut_ok = true;
			}
			else if (token == "max_computations") {
				if (nbc_ok) throw IOFileException("The number of maximum computations is defined more than once");
				token = tokeniser.GetNextToken(true);
				nb_comp_max = (unsigned int)Convert::StringToInt(token);
				nbc_ok = true;
			}
			else if ((token != "") || (!data.eof()))
				throw IOFileException("\"" + token + "\" is not a correct setup field");
		}
		cout << "=====================================" << endl;
		cout << "=             TSP Setup             =" << endl;
		cout << "=====================================" << endl;
		cout << "Number of bits per city    : " << nb_bpc << endl;
		cout << "Number of crossover points : " << nb_cross_pts << endl;
		cout << "Number of mutations        : " << nb_mutations << endl;
		cout << "Number of computations max : ";
		if (nb_comp_max==0) cout << "Infinite"; else cout << nb_comp_max;
		cout << endl;
		cout << "=====================================" << endl;
		// Close file
		data.close();
	}
	catch (exception &) {
		data.close();
		throw;
	}
}

void ClosePlugin () {
	// Output file stream
	ofstream data;
	// Save the solution to the associated file
	string file_name = program_name + ".sol";
	data.open(file_name.c_str(),ios_base::out|ios_base::app);
	// Test the file
	if (!data.is_open()) throw IOFileException("Cannot write to " + file_name);
	try {
		// Write data
		data << "Nb Cities     : " << distances.size() << endl;
		data << "Nb CO Pts     : " << nb_cross_pts << endl;
		data << "Nb Mutations  : " << nb_mutations << endl;
		data << "Best Solution : ";
		if (best_path.size() > 0) {
			data << "1 ";
			for (int i=0; i<(int)best_path.size()-1; ++i)
				data << (best_path[i]+2) << " ";
			data << (best_path[best_path.size()-1]+2);
		}
		data << endl;
		data << "Nb fitness computations : " << nb_computations << endl;
		data << "Best solution's length  : " << best_dist << endl;
		data << "Operators' history :" << endl;
		float ave;
		// Write operators fitness average
		for (unsigned int m=0; m<NB_OPS; ++m) {
			if (!mem[m].empty()) {
				ave = 0.f;
				for (unsigned int i=0; i<mem[m].size(); ++i)
					ave += mem[m][i].second;
				ave /= (float)mem[m].size();
				data << op_names[m] << " : Fit_ave = " << ave << endl;
			}
		}
		// End
		data << endl;
		// Close file
		data.close();
	}
	catch (exception &) {
		data.close();
		throw;
	}
	// Save the solution evolution data to the associated file
	data.clear();
	file_name = program_name + ".cvg";
	data.open(file_name.c_str(),ios_base::out|ios_base::app);
	// Test the file
	if (!data.is_open()) throw IOFileException("Cannot write to " + file_name);
	try {
		data << program_name << "_x = [ ";
		for (unsigned int i=0; i<record.size(); ++i)
			data << record[i].first << " ";
		data << "];" << endl;
		data << program_name << "_y = [ ";
		for (unsigned int i=0; i<record.size(); ++i)
			data << record[i].second << " ";
		data << "];" << endl;
		// End
		data << endl;
		// Close file
		data.close();
	}
	catch (exception &) {
		data.close();
		throw;
	}
	// Save the fitness results to the associated file
	data.clear();
	file_name = program_name + ".fit";
	data.open(file_name.c_str(),ios_base::out|ios_base::app);
	// Test the file
	if (!data.is_open()) throw IOFileException("Cannot write to " + file_name);
	try {
		for (unsigned int m=0; m<NB_OPS; ++m) {
			if (!mem[m].empty()) {
				data << op_names[m] << "_x = [ ";
				for (unsigned int i=0; i<mem[m].size(); ++i)
					data << mem[m][i].first << " ";
				data << "];" << endl;
				data << op_names[m] << "_y = [ ";
				for (unsigned int i=0; i<mem[m].size(); ++i)
					data << mem[m][i].second << " ";
				data << "];" << endl;
			}
		}
		data <<  endl;
		// Close file
		data.close();
	}
	catch (exception &) {
		data.close();
		throw;
	}
}

////////////////////////////////////////////////////////////////////////////////
// SC Callbacks
////////////////////////////////////////////////////////////////////////////////

// Initialise a path holder
bool Initialise (SCSystem & sys1, const SCSystem & context, SCSystem & sys2) {
	if (distances.size() == 0) throw AssertionException("No distance matrix defined for the TSP");
	// Create a list containing all the cities but the first,
	// sorted in enumeration order
	vector<unsigned int> list, path;
	for (int i=0; i<(int)distances.size()-1; ++i) list.push_back(i);
	// Create a random order
	vector<unsigned int>::iterator itor;
	while (list.size() > 0) {
		unsigned int idx = rand()%list.size();
		itor = list.begin(); itor += idx;
		path.push_back(*itor);
		list.erase(itor);
	}
	// Store it into the first schemata
	sys1.ReplaceSchema1(getStringPath(path,sys1.GetSchema1().size()));
	// Calculate the path's total distance and store it into the second schemata
	sys1.ReplaceSchema2(Convert::UIntToBinaryString(getTotalDist(path),sys1.GetSchema2().size()));
	// Send the solution holder into the computing space
	sys2.TakeFromSuperSystems(sys1, context);
	// Return OK
	return true;
}

// Keep best, PMX
bool KB_CO_PMX (SCSystem & sys1, const SCSystem & context, SCSystem & sys2) {
	return selectAndEvolve(sys1,context,sys2,KEEP_BEST,PMX);
}

// Keep best, MpOX
bool KB_CO_MOX (SCSystem & sys1, const SCSystem & context, SCSystem & sys2) {
	return selectAndEvolve(sys1,context,sys2,KEEP_BEST,MOX);
}

// Keep best, Swap Mutation
bool KB_CP_SWP (SCSystem & sys1, const SCSystem & context, SCSystem & sys2) {
	return selectAndEvolve(sys1,context,sys2,KEEP_BEST,SWP);
}

// Keep best, Move Mutation
bool KB_CP_MOV (SCSystem & sys1, const SCSystem & context, SCSystem & sys2) {
	return selectAndEvolve(sys1,context,sys2,KEEP_BEST,MOV);
}

// Keep best, Reverse Mutation
bool KB_CP_REV (SCSystem & sys1, const SCSystem & context, SCSystem & sys2) {
	return selectAndEvolve(sys1,context,sys2,KEEP_BEST,REV);
}

// Fitness proportional, PMX
bool FP_CO_PMX (SCSystem & sys1, const SCSystem & context, SCSystem & sys2) {
	return selectAndEvolve(sys1,context,sys2,FIT_PROP,PMX);
}

// Fitness proportional, MOX
bool FP_CO_MOX (SCSystem & sys1, const SCSystem & context, SCSystem & sys2) {
	return selectAndEvolve(sys1,context,sys2,FIT_PROP,MOX);
}

// Fitness proportional, Swap Mutation
bool FP_CP_SWP (SCSystem & sys1, const SCSystem & context, SCSystem & sys2) {
	return selectAndEvolve(sys1,context,sys2,FIT_PROP,SWP);
}

// Fitness proportional, Move Mutation
bool FP_CP_MOV (SCSystem & sys1, const SCSystem & context, SCSystem & sys2) {
	return selectAndEvolve(sys1,context,sys2,FIT_PROP,MOV);
}

// Fitness proportional, Reverse Mutation
bool FP_CP_REV (SCSystem & sys1, const SCSystem & context, SCSystem & sys2) {
	return selectAndEvolve(sys1,context,sys2,FIT_PROP,REV);
}

// Self-adaptation
//#define MIN(X,Y) (((X)>(Y))?(Y):(X))
bool AdaptOp (SCSystem & sys1, const SCSystem & context, SCSystem & sys2) {
	float	fit1 = getOpFit(sys1.GetKernel().substr(0,3)),
			fit2 = (&sys1 == &sys2)?fit1:getOpFit(sys2.GetKernel().substr(0,3));
	if ((fit1 < 0.f) || (fit2 < 0.f)) return false;
	
	/*
	float max = (fit1 + fit2) - MIN(fit1,fit2), wheel;
	wheel = ((float)rand()*max)/((float)RAND_MAX);
	if (wheel > fit1 - (MIN(fit1,fit2)/2.f)) {
	*/
	
	if (fit1 < fit2) {
		//if (fit1 <= 1.f) {
			sys1.ReplaceKernel(op_codes[op_wheel()],0,3);
			return true;
		//}
		//else return false;
	}
	else {
		//if (fit2 <= 1.f) {
			sys2.ReplaceKernel(op_codes[op_wheel()],0,3);
			return true;
		//}
		//else return false;
	}
}

////////////////////////////////////////////////////////////////////////////////
